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fNJ , A spherically symmetrical ultra-relativistic blast wave is not an attractor of a 

generic asymmetric explosion. Spherical symmetry is reached only by the time the blast 
wave slows down to non-relativistic velocities, when the Sedov-Taylor-von Neumann 
attractor solution sets in. We show however, that a directed relativistic explosion, 

CO ■ 

' with the explosion momentum close to the explosion energy, produces a blast wave 

with a universal intermediate asymptotic - a selfsimilar directed ultra-relativistic blast 
^Cj , wave. This universality might be of interest for the astrophysics of gamma-ray burst 

afterglows. 
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1. Relativistic Blast Waves and Gamma- Ray Burst Afterglows 

00 

o 

We give a selfsimilar solution describing the opening of a narrow ultra-relativistic blast wave. 
t^J- ■ This solution is an attractor (intermediate asymptotic) of a generic directed explosion. Here 

the "directed explosion" means an explosion with momentum P nearly equal to the explosion 
energy E (c = 1 here and below). Qualitatively, the directed blast wave solution has been already 
discussed by Rhoads (1999). 

•rH . 

A non-directed explosion (with E — P <^ E) does not have a universal intermediate asymptotic 
during the ultra-relativistic stage. The Blandford-McKee (1976) selfsimilar solution is an 
intermediate asymptotic only for a spherically symmetric explosion (Gruzinov 2000). A generic 
non-directed explosion acquires spherical symmetry only at the non-relativistic stage, when it 
asymptotes to the Sedov-Taylor-von Neumann attractor solution. 

Qualitatively, the non-universality of the Blandford-McKee selfsimilar solution follows from 
causality. The shock of the Blandford-McKee solution moves with Lorentz factor T oc i~ 3 / 2 . 
Consider a light signal propagating along the shock, starting at a polar angle 9 = 0. Then 
(dr/dt) 2 + (rd9/dt) 2 = 1, which, for T » 1, gives 9 = T . Two regions of the shock, separated 
by angle 9, do not talk to each other until the blast wave slows down to T ~ 6 . 

On the other hand, the selfsimilar solution derived here has a shock of opening angle 9 ~ 
The angular structure of the blast wave is not imposed by hand - it appears dynamically, from 
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an arbitrary initial state (satisfying the directed explosion requirement E — P <C E). Numerical 
simulations confirm the universality of the resulting selfsimilar solution. 

It is thought that ultra-relativistic blast waves are responsible for the observed gamma-ray 
burst afterglows (Piran 1999). Universality of the blast waves resulting from directed explosions 
should be of interest for astrophysics (Rhoads 1999). 

We give the main result in §2, derive it in §3, and describe the numerical simulations in §4. 



2. Selfsimilar Directed Ultra-Relativistic Blast Wave 

Selfsimilar directed blast wave is an axisymmetric flow with a shock. The shock is roughly 
parabolic. The head of the parabola moves at about the speed of light. The energy-containing 
part of the flow is a narrow shell behind the shock, with a growing opening angle. 

Consider a directed explosion of energy E which sends an ultra-relativistic blast wave into 
the medium of uniform density p. Let (r, z) be cylindrical coordinates, z is the axis of the directed 
explosion - meaning that the explosion occurs on the z-axis, with the explosion momentum in the 
positive z direction. 

Then the (proper) pressure p, the Lorentz factor squared q = 7 2 , and the radial (cylindrical 
radial ) velocity u behind the shock front are given by the following expressions 

p = e- 2t / T P(^ V )p, q = e- 2t ' T Q^ n ), u = e t / T U^r l ). (1) 



Here P(£,r]), Q(^,rj), U(£,r]) are dimensionless functions of dimensionless variables £, rj. 
Both the dimensionless variables and dimensionless functions are order unity quantities. 

The characteristic time r and the dimensionless variables are given by the following expressions 

where C is a dimensionless number of order unity. The solution applies for e~*/ T S> 1, that is for 
large negative t, when equation (1) gives an ultra-relativistic flow (with g> 1). 

We were able to calculate the dimensionless functions only approximately, with some ~ 10% 
accuracy. The flow fields P, Q, U, are given in figures (1, 2). Due to some inconvenient numerical 
features of the attractor solution (slow convergence of the energy integral, see §3.4), we were able 
to calculate only the approximate practial value of the dimensionless constant C - for the initial 
opening angle of the blast ~ 0.1, we find C « 1.5. 
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Fig. 1. — Left: Energy density E = iQP; E mSuX = 2.05, contours are E max over 2, 4, 8, 16. Right: 
Pressure P; P max = 0.83, contours are P max times 0.9, 0.7, 0.5, 0.3. 
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Fig. 2. — Left: Lorentz factor squared Q. ; Qmax = 0.63, contours are Q max times 0.9, 0.7, 0.5, 0.3. 
Right: Radial velocity U, contours are 1.0, 0.8, 0.6, 0.4, 0.2. 
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3. Basic Equations and Computations 

Once the good variables are identified (Blandford and McKee 1976), the computation is 
straightforward and very similar to the linearized case (Gruzinov 2000). In this section we 

1. write the relativistic hydrodynamics equations and the shock jump conditions in cylindrical 
coordinates, 

2. simplify the equations assuming the ultra-relativistic flow, 

3. find the right self-similar ansatz, and derive the equations describing the selfsimilar functions 

P(Z,v), Q(M, u(t,ri), 

4. perform the next similarity transformation, and find the asymptotic form of the functions 
P(£iW)i Q(£)?7)i U(£, 7]) in the far downstream region. 



3.1. Relativistic Hydrodynamics in Cylindrical Coordinates 

Relativistic hydrodynamics equations are written as d u T^ v = (Landau & Lifshitz 1987) 
. Here = (t, r) are Cartesian coordinates in the rest frame of the unshocked fluid. The 
energy- momentum is T^ u = {4ai^u v — g liU )p in the shocked fluid. The energy- momentum in the 
unshocked fluid, t^ u , has only one non-zero component, and without loss of generality in the 
final result, we take r 00 = 1. The shock position can be represented by the zero isosurface of 
a scalar field: <j) = at the shock. Then the boundary conditions at the shock follow from the 
hydrodynamics equations: 

T^dA = T^dvfa = 0. (3) 

The 4- velocity is = 7(1, v). In cylindrical coordinates, for a cylindrically symmetrical flow, 
v = vz + uf. Use the three-dimensional form of the hydrodynamics equations 

d T 00 + diT 0i = 0, (4) 
d T 0i + djT ij = 0, (5) 



and 



and project onto z and f. We get 

9 t [(4 7 2 - l)p] + 4d z [j 2 pv] + 4r~ 1 d r [r 7 2 H = 0, (6) 

4d t [j 2 pv] + d z [(4j 2 v 2 + l)p] + Ar^drlr^pvu] = 0, (7) 

4dt h 2 pu] + 4d z [-f 2 pvu] + 4r" 1 <9 r [r^pu 2 ] + d r p = 0. (8) 



- 5 - 



Parameterize the shock front as (f> = z — z s (t, r) = 0. Then eq.(3) gives: 

4j 2 p[d t z s -v + ud r z s ] - pd t z s = d t z s 

<i^ 2 v[dtz s — v + ud r z s ] — 1 = 0. 
4j 2 u[dtz s — v + ud r z s ] + d r z s = 0. 



(9) 
(10) 
(11) 



3.2. Ultra-Relativistic Hydrodynamics in Cylindrical Coordinates 

Introduce a coordinate grid moving at the speed of light: x = t — z (giving the following 
replacements: dt — * dt + d x , d z — > —d x ). Parameterize the shock position as z s = t — x s (t,r). 



Denote j 2 = q. Then, 



> - l _iV_( 1 +^ +Q(f 3 ) , 



2q 



8q 2 



(12) 



Here and in what follows we assume g > 1 and qu 2 = 0(1). The first condition, q 3> 1, means 
that the flow is ultra-relativistic. The scaling qu 2 = 0(1) is confirmed by the result. 



Keeping the two leading orders in q 1 , we get 



dt[(4q - l)p] + d x 
d t [(4q-2-2qu 2 )p]+d x 



,2\2 



l + 2qu 2 + ^±^ ]p 



l + 2gw — 

2g 



2q 

2\2\ 



4r 1 3 r [rg'pu] = 0, 



P 



Ad t [qpu] + 2<? x [(l + qu 2 )pu] + 4r _1 d r [rgpu 2 ] + d r p = 0. 
In the leading order in q -1 , these can be written in the following final form 
dt(qp) + jd x [(l + 2qu 2 )p] + r^drirqpu] = 0, 
8 t [(l + 2<?u 2 )p] + d x [(l + qvffq^p] + 2r~ 1 a r [r(l + gn 2 )pn] = 0, 



d t u + 



1 + gti 



1 



d x u + - — d x p + - — d r p = 0. 
Iq 4qp 4qp 



The boundary conditions, at x = x s : 

Q 



1 _4 
4^ + 2(8^)2' P "3 9 ' 



(13) 



+ r _1 a r [r(4^ - 2 - 2gu 2 )H = 0, (14) 



(15) 

(16) 

(17) 
(18) 

(19) 
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3.3. Selfsimilar Solution 

The ultra-relativistic blast wave equations (16)- (19) admit the following scalings: 

r<xt 1/2 x 1/2 , pozqoctx' 1 , it oc t" ! / 2 x 1/2 (20) 

In D spatial dimensions, the energy of the flow behind the shock corresponding to these scalings is 

E ~ xr u L pq on t 2 x 2 . (21) 

During the ultra-relativistic stage the energy of the flow should be close to the explosion energy, 
meaning that the energy (21) should be time- independent. This enforces the following scaling for 
the directed blast wave in D spatial dimensions: 

3+-D 

xcxt~. (22) 
For D = 3 spatial dimensions this algebraic dependence degenerates into the exponential. 
We therefore want to find the following selfsimilar solution of the system (16)- (19): 

p = e- 2t P(Z,7 1 ), q = e- 2t Q(d,ri), u = e t U(^, V ), (23) 

with 

£ = e ~ 2t x, 7] = e~V. (24) 
Using the ansatz (23), (24) in the equations (16)-(19), one gets 

(_ 4 _ - V d v )(QP) + ^[(1 + 2QU 2 )P] + ri-%[riQPU] = 0, (25) 

(-2 - 2£d c - 7]d v )[(l + 2QU 2 )P] + <%[(1 + QU 2 ) 2 Q- l P] + 2^ 8^(1 + QU 2 )PU] = 0, (26) 

(1 - 2^ - + 

The boundary conditions at £ = £ s (v) : 



(i - m - ,a„)u + l Agfw + + = »■ (^) 



q = P = ^ ' " = e " (28> 

The explosion energy is 

E = 8ir J rjdrj QP, (29) 

— 1/3 

giving the dimensionless constant C of §2: C = -E 

Equations (25)-(28) still have a one-parameter scaling group: 

£^A 2 £, r?^Ar/, Q -> X~ 2 Q, P -» A~ 2 P, U ^ XU. (30) 

"Solution" shown in figures 1, 2 was normalized by the condition £ s (0) = 0.2. This normalization 
gives order unity values for all other quantities. 

The right way to find the attractor is to solve the equations (25)-(28). This is not what we 
have done in this paper. Our "solution" was obtained by a direct numerical simulation of the 
axisymmetric relativistic hydrodynamics equations (§4). 
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3.4. Selfsimilar Selfsimilar Solution 

Far from the head of the blast wave, at large values of n, the selfsimilar attractor equations 
have a selfsimilar solution corresponding to the power law shock position 

6 = *A (31) 

If (3 < 2, the energy integral (29) diverges. This does not mean that we may automatically 
dismiss such solution. It might rather indicate that the attractor (1) is never fully filled up. We 
therefore looked for such solutions. It was found that only the (3 = 2 solution (with logarithmically 
divergent energy) exists. 

For (3 > 2, the scaling (31) enforces the following scalings of the fields: 

Q = vf 2 ^-^q{x), P = V - 2(/3 -Vp(x), U = r} (} - 1 u{x), x = ^, (32) 
with the boundary conditions on the shock: 

q{l) = w p{l) = w u{l)=l3 (33) 

The functions q, p, u are obtained from the system of ordinary differential equations: 

-[(1 + 2qu 2 )p}' - (3/3 - 4)qpu - (3x{qpu)' = (34) 

[(1 + <?ti 2 ) W - (2/3 - 4)(1 + qu 2 )pu - 2(3x[(l + qu 2 )pu}' = (35) 

2{l + qu 2 )u' - 2{(3 - I) + (u - (3x)- =0 (36) 

p 

Numerical integration of equations (34)- (36) then shows that solutions exist only for 
(3 < 2.032. But we were unable to determine the true value of (3 theoretically. The direct numerical 
simulations (§4) give marginal evidence that j3 > 2 (the case (3 = 2 has positive u for all x, the 
case (3 > 2 has negative u for large enough x, we do see negative velocities in the far downstream 
region). We therefore tentatively conclude that 2 < [3 < 2.032. 

Since (3 is so close to the minimal value (3 = 2, the energy integral (29) converges very slowly. 
As a result, only extremely directed explosions can reach the pure attractor stage. Explosions 
with initial opening angles of order 0.1 or 0.01 produce only a partial attractor - the attractor is 
being filled up throughout the entire quasi-selfsimilar stage. This should give small (algebraic) 
corrections to the exponential decay laws of maximal energy and pressure of the pure solution (1). 
Direct numerical simulations (§4) give (3 « 2, to about 10% . 
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4. Numerical Simulation 

Simple Lax scheme with 800x800 resolution was used to simulate the axisymmetric relativistic 
hydrodynamics in cylindrical coordinates eq.(6-8). Coordinate mesh was moving at the speed of 
light - this allows to simulate just a small region around the energy containing part of the flow. 

Different initial configurations were tried and seen to produce approximately equal attractors. 
The solution shown in fig. (1,2) used the following initial condition. We started with a blob of 
energy moving in the z direction with the maximal Lorentz factor equal to 20. Maximal initial 
pressure in the blob was 400; the blob was moving into the medium of energy density equal to 1. 
Maximal initial energy density in the blob was 4pq = 640, 000. The simulation was stopped when 
the maximal energy density of the flow dropped to about 2500. Initially, the blob had the r to z 
size ratio equal to 1:20. 

Figures (1,2) show the final state of the blob, in the appropriately rescaled variables. The 
dimensionless constant C was calculated from the initial energy of the blob and the inferred time 
constant r measured from the numerical solution. 

It was found that the blob does approximately evolve toward the attractor solution given by 
(1), (2). Namely (to about 10% accuracy): 

1. Maximal energy density e max of the flow decreased just somewhat faster than exponentially. 

2. Position of the shock x s grew almost exponentially. 

3. e m ax and x s were approximately related by e max oc x~ 2 in agreement with (23), (24). 

4. Total energy per logarithmic interval of energy density was approximately constant for low 
energy densities, corresponding to [3 « 2. 

5. The time constant r inferred from the rate of change of maximal energy density e max , was 
approximately equal to r = 8x s <7 max , in agreement with (28). 

We should not have expected a better than 10% agreement. First, Lax method does not give 
accurate maximal values (which theoretically occur at the shock). At resolution of 400x400, some 
of the inferred values change by about 10% . Second, we are simply not ultra-relativistic enough. 
We stop at the maximal Lorentz factor of about 5. A fair fraction of the flow has too small Lorentz 
factors. Finally, because of the slow convergence of the energy integral (§3), the attractor is being 
filled up all the time - hence faster than exponential decrease of the maximal energy density. 

I thank Andrew MacFadyen for discussions. This work was supported by the Davide and 
Lucile Packard Foundation. 
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